function [t, Y] = rk4(fun, tspan, y0, h)
    % 使用经典的四阶 Runge-Kutta 方法求解常微分方程
    % 输入：
    %   fun: 函数句柄，形式 dydt = odefun(t,y)
    %   tspan: [t0, tf]
    %   y0: 初始条件向量
    %   h: 步长
    % 输出：
    %   t: 时间向量
    %   Y: 状态矩阵，每行对应一个时间点的解
    
    t0 = tspan(1);
    tf = tspan(2);
    t = t0:h:tf;
    if t(end) < tf
        t = [t, tf];
    end
    n = length(t);
    Y = zeros(n, length(y0));
    Y(1,:) = y0(:)';
    
    for i = 1:n-1
        ti = t(i);
        yi = Y(i,:)';
        h_step = t(i+1) - t(i);  % 考虑最后一步可能小于 h
        k1 = fun(ti, yi);
        k2 = fun(ti + h_step/2, yi + (h_step/2)*k1);
        k3 = fun(ti + h_step/2, yi + (h_step/2)*k2);
        k4 = fun(ti + h_step, yi + h_step*k3);
        Y(i+1,:) = yi + (h_step/6)*(k1 + 2*k2 + 2*k3 + k4);
    end
end